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In this paper a new procedure to smooth out the initial energy densities of hydrodynamics is 
employed to show that the initial spatial eccentricities £ m , n , which drive the final state flow har¬ 
monics v n , are remarkably robust with respect to variations of the underlying scale of initial energy 
density spatial gradients, A, in nucleus-nucleus collisions. For = 2.76 TeV Pb+Pb collisions, 
the £m,U s (across centrality classes) change by less than 10% if the scale of fluctuations is varied 
from 0.1 to 1 fm. We show, using the 2+1 Lagrangian hydrodynamic code, v-USPhydro, that this 
robustness is transferred to the final r> n ’s computed within event by event viscous hydrodynamics. 

This indicates that the flow harmonics in nucleus-nucleus collisions are not yet particularly sensitive 
to the underlying microscopic sub-nucleon physics below the confinement scale. On the other hand, 
the eccentricities of top 1% high multiplicity y/s = 5.02 TeV p+Pb collisions are found to be very 
sensitive to sub-nucleonic scale fluctuations, which should be contrasted with the robustness found 
in peripheral Pb+Pb collisions with the same multiplicity. 
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I. INTRODUCTION: A TALE OF SCALES IN NUCLEUS-NUCLEUS COLLISIONS 

After more than a decade of intense investigation there is by now large experimental evidence that the local initial 
energy density of the matter formed in ultrarelativistic heavy ion collisions is sufficiently large enough to produce a 
deconfined plasma of quarks and gluons The main experimental signatures of this new state of matter can be 

divided into two groups. On one hand, the quark-gluon plasma is color opaque with respect to the propagation of 
partonic jets, which leads to the experimentally observed suppression of high pp single inclusive hadrons known as jet 
quenching [MS] (see 0 for a recent theoretical discussion). On the other hand, the large anisotropic flow coefficients 
of low transverse momentum hadrons, pr < 2 GeV, indicate that there is a significant degree of collectivity in the 
system, which is generally consistent with relativistic hydrodynamic calculations in the “perfect fluid” limit where 
viscous effects are minimal (for a recent review see QZIj). This suggests that the non-Abelian quantum fields present 
in the initial stage of the collisions may behave incoherently at length scales (at least) on the order of the size of 
a large nucleus in a way that is compatible with the strong assumptions behind relativistic hydrodynamic behavior 

PH EH- 

However, just as the initial inhomogeneities in the early universe led to the formation of large scale structures after 
expansion [13], the hot and dense matter created in heavy ion collisions displays energy density fluctuations of different 
characteristic sizes and physical origins. For instance, in the context of heavy ion collisions one may consider the size 
of a large lead nucleus Rp 5 ~ 7 —10 fm as a macroscopic nuclear length, the size of a nucleon ~ 1 fm as an intermediate 
mesoscopic scale, while smaller length scales of the size of the inverse of the saturation scale ^ 1 /Q S < 0.1 fm can be 
considered as the microscopic sub-nucleon regime. The macroscopic nuclear scale is essentially “geometrical” in the 
sense that it can be varied by changing the type of nucleus, the final hadron multiplicity yield, and the center-of-mass 
collision energy per nucleon. The mesoscopic length scale is of the order of the inverse nonperturbative 1/A qcd scale 
and at that scale one cannot yet resolve the color field configurations inside the nucleons. Finally, the microscopic 
sub-nucleon scale resolves the internal color structure of the nucleon and some of its features can be understood at 
weak coupling (though in a nonperturbative manner, for a review see [14]). Therefore, ultrarelativistic collisions of 
heavy nuclei involve QCD phenomena associated with characteristic length scales that can vary by approximately 
three orders of magnitude. At the moment, there is no effective theory that can describe phenomena at these different 
scales in a single, consistent framework. 

While variations in the macroscopic nuclear scale have been studied at length both theoretically and experimentally 
over the years m , the interplay between mesoscopic and microscopic sub-nucleon scales and their effects in the 
calculation of heavy ion collision observables, such as the anisotropic flow coefficients, still require further investigation. 
Given the success of relativistic hydrodynamics to describe anisotropic flow data in nucleus-nucleus collisions, it is 
reasonable to assume that hydrodynamics may be an adequate effective theory framework to be used in such a study. 

Current numerical simulations of relativistic hydrodynamic employ Israel-Stewart-like models where the shear 
viscous contribution to the energy-momentum tensor, 7 ^, obeys an independent differential equation of the relaxation 
typ<0 This relaxation process towards the universal Navier-Stokes regime is characterized by the shear relaxation 
time transport coefficient r n . For dilute gases described by the relativistic Boltzmann equation, the 14-moments 
approximation for a massless, single component gas (assuming classical statistics) with constant cross section gives 

01 EH 

5 p 

T?r = TF 

where T is the temperature, 77 is the shear viscosity, and s is the local entropy density. For instance, if we assume 
that the largest temperature achieved in central Pb+Pb collisions at the LHC is, say T ~ 400 MeV, for p/s = 0.1 the 
formula above gives r n ~ 0.25 fm. An important aspect of this formula is that even for a constant p/s the relaxation 
time may vary significantly in space and time according to the local temperature profile. 

However, as discussed in m, in Israel-Stewart (IS) models is the relevant microscopic scale associated with 
shear stress and, as such, it provides a lower bound on the scales at which one expects this hydro dynamical model 
to provide accurate result^] In fact, microscopic sub-nucleon processes that involve scales shorter than are not 
properly taken into account in IS models and their description requires the introduction of other degrees of freedom 
beyond those already present in IS theory. This has been shown to be the case in [18], in a systematic manner, 
using the relativistic Boltzmann equation as the underlying microscopic theor}0 Alternatively, the current success 


1 This implies that in these models one needs to provide not only initial conditions for the standard hydrodynamic fields, such as the 
energy density e (or pressure P) and flow velocity u^, but also must be independently known at the initial time to. Note that in 
this paper we do not consider effects from other conserved charges such as baryon number. 

2 Clearly, if the bulk viscosity contribution to the energy-momentum tensor, II, is included in the dynamics one has to take into account 
its corresponding relaxation time rn as well. 

3 At strong coupling the equations of motion that describe the evolution of the shear stress tensor are not of relaxation-type [17111191120] 

and, in this case, IS theory fails to properly describe transient fluid dynamical phenomena. 
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of anisotropic hydrodynamics in describing strongly dissipative processes with large momentum space anisotropy 
requires the introduction of new degrees of freedom in the dynamics that are not the usual hydrodynamic fields (for a 
recent study, and references therein, see [21]). Therefore, IS-like theories do not correctly resolve scales below and 
this intrinsic limitation of this effective theory must be taken into account in numerical simulations. In fact, such a 
limitation becomes especially relevant in event by event hydrodynamic simulations in which the hydrodynamic fields 
may possess large spatial gradients at early times. 

The presence of large spatial gradients can be verified by computing the Knudsen number, iT n , which is roughly 
speaking the product of the relevant microscopic scale (in our case, r w ) and a spatial gradient of the hydrodynamic 
fields (in this sense, it is clear that K n should be small for hydrodynamics to be valid). Using the gradients of e and 
uv one can construct a handful of different Knudsen numbers (or fields since these quantities vary in space and time). 
In [22 the regions where K n = 0.5 were considered to define the edge of the validity of the hydrodynamic description 
and it was found that the Knudsen number K u q = r n 0 associated with the local expansion rate, 0 = and 

the Knudsen number K ne = r^y/V ^e^V^e/e, associated with the gradient of the energy density, were particularly 
relevant to assess the overall validity of the hydrodynamic description of the QGP formed in ultrarelativistic collisions. 

The different length scales discussed above and the assumed region of applicability of IS hydrodynamics are sum¬ 
marized in Fig. [l] where A is the smoothing scale of the hydrodynamic simulations. We have also included snapshots 
of a typical initial energy density profile that show how the same event changes as the smoothing scale increases (this 
will be discussed in detail in the next section). We also note that hydrodynamics cannot account for phenomena 
that occur at microscopic sub-nucleon scales at which the effects from coherent, color fields should not be neglectecQ 
Clearly, as the smoothing scale A increases, the local gradients ~ 1/A decrease (and so do the Knudsen numbers), and 
only the meso and macroscopic nuclear scale regimes are accessible in the hydrodynamic description. 

The smoothing scale provides a natural way to characterize the different models for the initial conditions of hy¬ 
drodynamics. For instance, models where the sources of fluctuations are related to the positions of the uncorrelated 
nucleons in the colliding nuclei, such as the Monte Carlo (MC) Glauber model [331135] and the MCKLN model [36] . 
cannot properly resolve microscopic sub-nucleon scales (in the sense discussed above). On the other hand, the IP- 
Glasma model [37, or the Dumitru-Nara BK model [35] depend on the physics at microscopic length scales of 0(1/Q s ), 
which in turn produces large energy density fluctuations and structure even at very small scales. It is not clear at the 
moment how these very short wavelength fluctuations affect observable quantities related to hydro dynamical behavior 
such as the anisotropic flow coefficients. 

In this paper we study the dependence of anisotropic flow coefficients on variations in the smoothing scale present 
in the initial conditions used in hydrodynamic simulations. The smoothing procedure is based on a Lagrangian 
description of hydrodynamics but it can also be employed to smooth out initial conditions that may later be used 
in (Eulerian) grid codes as well. The method can be used to smooth out initial conditions from both Glauber-like 
models as well as saturation based models such as IP-Glasma. In our approach, after going through this “UV filter” 
defined by the smoothing scale A, the initial condition only contains scales of size A and above: shorter scales are 
consistently smeared out. After this procedure, one can perform event by event viscous hydrodynamic simulations to 
investigate how smoothing out the short scales in the initial conditions may affect the hydrodynamic response of the 
fluid to the initial spatial inhomogeneities. 

We consider here MCKLN initial conditions [36] for y/s = 2.76 TeV Pb+Pb collisions at the LHC and investigate 
how the anisotropic flow coefficients computed using the 2+1 viscous hydrodynamics code v-USPhydro f39h4T] are 
affected by changes in the smoothing scale. While the effects of different length scales present in the initial state 
have been investigated using different approaches |42h5T] , an important novelty of our work is that the initial spatial 
eccentricities for nucleus-nucleus collisions (across centrality classes) do not change appreciably when the smoothing 
scale is varied in the mesoscopic regime ranging from 0.1 fm to 1 fm (a similar result was found in the context of 
the PHSD transport code [51]). This result illustrates that, at least for the case of MCKLN initial conditions (and 
MCGlauber events), length scales below the mesoscopic scale (^ 1 fm) given by the uncorrelated fluctuations in the 
position of the nucleons have a very small effect in the final azimuthal anisotropies in nucleus-nucleus collisions. Our 
results qualitatively agree with the conclusion of Ref. [46] (obtained within mode-by-mode hydrodynamics) that the 
azimuthal anisotropies are much more sensitive to intermediate length scales (defined in the mesoscopic regime in our 
notation) than to scales deep in the microscopic sub-nucleon scale regime. Once the smoothing parameter enters the 
macroscopic nuclear scale regime, the eccentricities begin to change and if the process is continued even further all 
the structure is lost and nontrivial azimuthal anisotropies, such as triangular flow, would not be observed. 

We also study how the variations in the smoothing scale affect the initial spatial eccentricities of top 1% high 
multiplicity y/s = 5.02 TeV p+Pb collisions computed using the recently developed Trento code [52]. We find that 


4 We note that this transition from a color field dominated regime to hydrodynamics at large scales has also been extensively studied in 
the past in the context of jets E312]. In fact, even in the case of a static and uniform medium, the small region near a fast moving 
colored parton is too influenced by the local color fields to display hydrodynamic behavior. 
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FIG. 1. (Color online) Different length scales involved in the ultrarelativistic collisions of heavy nuclei separated according 
to the smoothing scale A. The energy density snapshots on the right illustrate how the initial condition changes with the 
smoothing scale. 


the eccentricities in p+Pb collisions are particularly sensitive to variations in the smoothing scale in the sub-nucleon 
regime. This sensitivity to the fluctuation scale in the mesoscopic regime is not seen in 65 — 70% peripheral yfs = 2.76 
TeV Pb+Pb collisions even though their multiplicity is comparable to high multiplicity p+Pb collisions [53] , 

Also, in this paper we study how variations in the smoothing scale on an event by event basis affect K u q and the 
inverse Reynolds number, Re~ x = W jp [l|| in Pb+Pb collisions. While Re~ x remains generally small for most 

of the evolutiorj^] in the mesoscopic regime, K u q is significantly enhanced by the spatial inhomogeneities present in 
event by event simulations. However, given the robustness of the azimuthal flow coefficients with respect to variations 
of A in the mesoscopic regime in nucleus-nucleus collisions, hydrodynamic events with large local Knudsen numbers 
can be formally “evolved” in A to reduce their local Knudsen number while keeping final hydrodynamic observables 
nearly unchanged. 

This paper is organized as follows. In Section [IT] we explain our smoothing procedure and show how it affects 
the initial energy density and spatial eccentricities in Pb+Pb collisions. In Section m we give the details about 
our hydrodynamic simulations and show how variations in the smoothing parameter affect the Knudsen and inverse 
Reynolds number in an event by event simulation. We show in Section [IV] our results for the anisotropic flow coefficients 
in Pb+Pb collisions and make a comparison to LHC data. Our study on the dependence of the eccentricities with 


5 Our trivial choice for the initial conditions of the shear stress tensor at the initial time to, i.e, 7t^(to, r) = 0 plays a big role to ensure 
this condition. 
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the smoothing scale in p+Pb collisions is done in 
and outlook. 

Definitions : We use a mostly minus metric 
where 


Section [V] We finish the paper in Section [VI] with 


= (1, —1, —1, —r 2 ) expressed in Milne coordinates 


our conclusions 
= (r,x,y,<;) 



~ ^ , 
t -\- z 
t — z 


( 2 ) 

(3) 


are the propertime and spacetime rapidity (expressed in terms of the standard Minkowski coordinates), respectively. 
Throughout this paper we use natural units where h = c = ks = 1. Coordinates in the transverse plane are denoted 
as r = (x, y). 


II. SMOOTHING OUT INITIAL CONDITIONS FOR HYDRODYNAMICS 

Different assumptions regarding the underlying physics of the initial conditions for hydrodynamics in heavy ion 
collisions can lead to very different initial energy density profiles. For instance, MC Glauber-based models [34: produce 
much smoother (though of course still largely inhomogeneous) event by event energy density profiles than those found 
in IP-Glasma m- It is of course of interest to check how hydrodynamic observables (such as the anisotropic flow) 
change with the spatial resolution present in the initial conditions. 


A. Definition of the smoothing procedure 

Initial condition models give (in one way or the other) the expectation value of the energy-momentum tensor of the 
matter, T Mzy (ro, r), at an initial time To for each evenf[^] The idea here is to use some type of UV filter to systematically 
remove (on an event by event basis) from the initial T^ u fluctuations of wavelength below a given smoothing scale A. 
We shall denote the energy-momentum tensor that has gone through this A—filter at the initial time by T^ v (to, r; A), 
i.e., 

(to , r) - (to , r; A). (4) 

To study how variations in the smoothing scale change the initial energy-momentum tensor one needs to solve a 
smoothing flow functional equation 

r; A) = T U u [T a p(T 0 , r; A); A], (5) 

where is a functional of the energy-momentum tensor. Clearly, the specific form of T^ v depends on the details 
regarding the implementation of the UV-filter. While the derivative in (J5| is taken only with respect to the smoothing 
parameter, we expect that with the smoothing flow the spatial structures in the initial condition are smoothed out 
(or diffused) up to the minimal scale A and only the large scale properties of the energy-momentum tensor surviv^J 
Thus, given an initial T^ v that displays structure in all the regimes displayed in Fig. |T} the smoothing ffovj^] procedure 
explained above can be used to systematically remove microscopic sub-nucleon scales and bring the system towards 
the meso or macroscopic nuclear scale regimes. 

A simple UV-filter that works as described above is the one already currently used in Lagrangian methods to 
solve the hydrodynamic equations, such as the Smoothed Particle Hydrodynamics (SPH) formalism [MI [56], which 
has been successfully employed for more than a decade to study event by event simulations of numerical relativistic 
hydrodynamics in the context of heavy ion collisions within NeXSPheRIO [571465] and also more recently in the v- 
USPhydro code [39H4T] . The SPH implementation used in v-USPhydro was discussed in detail in [39! and we shall 
review below only the points needed to define the UV-filter discussed above. 


6 In this paper we only discuss longitudinally boost invariant systems so no dependence on the spacetime rapidity is included. Clearly, 
our discussion can be generalized to include a nontrivial rapidity dependence as well, though in this case one may want to distinguish 
the smoothing scale in the transverse plane from that used in the rapidity direction. 

7 One would be remiss to not point out here a loose resemblance with Ricci flow In the latter, the metric of a manifold is evolved 
in a process formally analogous to heat diffusion to smooth out irregularities in a way that the knowledge of the geometric structure of 
the manifold can lead to useful global topological information. In our case involving hydrodynamics, the energy-momentum tensor is 
evolved in A, removing small scale irregularities, while the global overall structure is preserved. 

8 We would like to stress that the “flow” in the smoothing procedure explained here has to do with variations in the underlying scale of 
the initial condition - no hydrodynamic evolution is done at that point and, thus, one should not identify this smoothing flow with the 
truly dynamical hydrodynamic flow. 
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For simplicity, let us assume that the initial energy-momentum tensor is diagonal (i.e., there is no initial hydrody¬ 
namic flow or viscous contributions) such that we only need to consider the local energy density and pressure, which 
are assumed to be related via the equation of state (EOS). Of course, given the energy density the usual thermody¬ 
namical identities allow one to compute the local entropy density s(to, r) and temperature T(to, r) profiles. The main 
idea is that the initial energy density can always be reconstructed as follows: 

e( r o, r ; A) = ^ e a (ro) W ^A a), (6 ) 

rv = 1 ' ' 


where e a (ro) denotes a discrete set of variables associated with different points in the transverse plane, {r a , a = 1, ..., N] 
where one places a piecewise continuous distribution function W (|r — r a |/A; A). The kernel W has finite support^] 
given by A and, thus, it strictly vanishes for |r|/A sufficiently larger than unit. The variables e a (tq) must be carefully 
chosen in order to better describe the structures in the initial condition for a given value of y] In this work we use 
the cubic spline kernel employed in [39l [4Tj given by 


where 
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(7) 


(8) 


One can see in Fig. [2] that / strictly vanishes when |r|/A > 2. Also, from © one finds that (after filtering) the 
local gradient de ~ 0[ 1/A), which implies that as A increases the energy density becomes smoother. However, it is 
important to notice that this filter somewhat preserves the relative peaks and valleys in the energy density. In fact, 
let r a be the position of a peak in the energy density. One can see that an increase in A does not change the location 
of the peak (e a and r a remain the same) though its magnitude does decrease. In fact, for the implementation of the 
smoothing flow used in this paper one can see that 


d 

dX 


e(r 0 ,r; A) 


2 

A 


e(r 0 ,r; A) 



(9) 


9 The kernel is a delta sequence, i.e., lim^o VF (|r|/A; A) = 5(r) and is also normalized, i.e., f d 2 rW (|r|/A; A) = 1. 
10 See the discussion in m based on the variational principle. 
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a.) A = 0.1 fm b.) A = 0.3 fm 



FIG. 3. (Color online) Effects of the smoothing procedure on the energy density profile of an MCKLN event in the 0 — 5% 
centrality class of = 2.76 TeV Pb+Pb collisions at the LHC. The following smoothing parameters were used: a.) A = 0.1 
fm, b.) A — 0.3 fm, c.) A = 1 fm, and d.) A = 3 fm. 


which illustrates the points discussed above. Also, even though the first term on the right hand side of the equation 
above (—2e/A) will be the same for any normalized kernel, the relative importance of the second term with respect 
to the first may depend on the properties of the kernel itself. 

We remark that even though the discretization procedure used here was originally devised to be used in a Lagrangian 
approach to hydrodynamic^] the procedure explained in this section can be employed to smooth out initial energy 
density profiles that can be later used as input in Eulerian codes as well. 


B. Smoothing out MCKLN events 

Here we apply the procedure discussed above to smooth out initial energy densities generated by the MCKLN 
code of [36] for yfs = 2.76 TeV Pb+Pb collisions at the LHC. In Fig. [ 3 ] we show how the energy density of an event 
belonging to the 0 — 5% centrality class changes if one increases the smoothing parameter from A = 0.1 to 3 fm. One 
can see that small inhomogeneities are progressively smeared out with increasing A until only the large scale energy 
density structure in panel d.), corresponding to the macroscopic nuclear scale regime mentioned in Fig. [l] is visible. 


11 In fact, r a can be readily interpreted as the Lagrangian coordinates in a subsequent hydrodynamic evolution as done, for instance, in 

E3HD- 



1 . Initial eccentricities 


In this context, it is instructive to check how the initial spatial eccentricities vary with A. The eccentricities 
computed event by event, £ m?n , acquire a A dependence due to their definition in terms of the initial energy density 


iW 1 


.(A) _ 


f d 2 r r m e m< ^ e(ro, r; A) 
f d 2 r r m e(r 0 , r; A) 


( 10 ) 


where r = |rL 0 = tan 1 (t//x), and ^ m?n are the initial angles. The dependence of £2 ,n {n = 2 ,... , 6) with A can be 
seen in Fig. El for two different centrality classes (0 — 5% on the left and 20 — 30% on the right FM In this plot we 

show £2 , n = y (^2 ,n) with the average being performed over 150 events in each centrality class. One can see that the 

eccentricities only change appreciably after one enters the macroscopic nuclear scale regime A > 1 fm - this is the case 
for both centrality classes (we have also checked that this is the case for 0 — 1% and 65 — 70%). This indicates that 
the scale that sets the upper limit of the mesoscopic scale is indeed of the order of 1/A qcd ~ 1 fm. However, note 
that higher order eccentricities, such as £ 2 , 6 , begin to change for values of A smaller than £ 2,2 does. We remark that 
the robustness of the eccentricities with respect to variations of A found here are consistent with the previous general 
analysis performed in [66] . 

The relative variation of the eccentricities A£ m?n = 100 [£ m , n (Ao) — £™,n(A)] /£ m ,n(A 0 ) is displayed in Fig. [6} which 
shows that the lowest order eccentricities indeed change by less than 10% when A varies from 0.1-1 fm (the higher 
order eccentricities are indeed more sensitive to smoothing, as expected). In summary, for MCKLN initial conditions, 
only when the smoothing scale is larger than 1 fm the essential spatial structures that contribute significantly to 
the eccentricities (and to the final azimuthal anisotropies) are washed out and the eccentricities themselves start to 
decrease significantly. 



A [fm] A [fm] 


FIG. 4. (Color online) Dependence of the initial eccentricities £ 2 ,n = y (^ 2 ,n) (f° r n = 2 — 6) on the smoothing parameter, A, 

for = 2.76 TeV Pb+Pb MCKLN events. The left plot was computed using events in the 0 — 5% centrality class while on the 
right events in the 20 — 30% centrality class were used (the average was performed over 150 events for each centrality class). 


The same conclusion can drawn about the other eccentricities shown in Fig. [5j which were previously shown to be 
relevant in the prediction of the final flow harmonics (see for instance, EOT). Since the fluctuations in MCKLN 
have the same origin as in MCGlauber models (i.e., the random positions of the nucleons), we believe that the results 
found here concerning the dependence of eccentricities with the smoothing scale should also hold for MCGlauber initial 
conditions. It would be interesting to investigate how the eccentricities computed using the much finer IP-Glasma 
initial conditions m vary with A. One expects that, because of the presence of sub-nucleonic fluctuations in this 
type of model, its eccentricities should be more sensitive to the smoothing scale and significant changes may be found 
already for smaller values of A than those found here in MCKLN (i.e., the upper limit of the mesoscopic regime occurs 
at A < 1 fm). 
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The initial energy densities before smoothing are defined on a square lattice of size Aq = 0.06 fm. 
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FIG. 5. (Color online) Dependence of the other initial eccentricities £ m , n = yj(im,n) on the smoothing parameter, A, for 
yfs — 2.76 TeV Pb+Pb MCKLN events. The left plot was computed using events in the 0 — 5% centrality class while on the 
right events in the 20 — 30% centrality class were used (the average was performed over 150 events for each centrality class). 
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FIG. 6. (Color online) Relative variation of initial eccentricities As 2 , n = 100[£2,n(Ao) — £ 2 ,n(A)] /e 2 ,n(Ao) on the smoothing 
parameter, A, for yfs = 2.76 TeV Pb+Pb MCKLN events. The left plot was computed using events in the 0 — 5% centrality 
class while on the right events in the 20 — 30% centrality class were used (the average was performed over 150 events for each 
centrality class). 


III. DETAILS ABOUT THE HYDRODYNAMIC EVOLUTION AND RESULTS FOR THE KNUDSEN 
AND INVERSE REYNOLDS NUMBERS IN NUCLEUS-NUCLEUS COLLISIONS 

We used the v-USPhydro code 391441] to perform hydrodynamic simulations using as initial conditions the y/s = 2.76 
TeV Pb+Pb MCKLN events in the 0 — 5% and 20 — 30% centrality classes discussed in the last section. The current 
version of v-USPhydro solves the energy-momentum conservation equations in Milne coordinates 

= 0 =*► f (tTH + = 0, (11) 

where 

= \a va (d^gaf) + dpg a ^ - d a g^) ( 12 ) 

is the Christoffel symbol computed using the Milne metric g^ v = (1, —1, —1, —r 2 ). Though bulk viscosity has been 
shown to be relevant in event by event hydrodynamic simulations urn EH EHHZD], for simplicity in this paper we 
focus solely on shear viscous effects and set the bulk viscosity to zero. Contributions from other 2nd order hydrody¬ 
namic transport coefficients besides r nj such as those computed either in kinetic theory [18] or in strongly coupled 
nonconformal holography ED, are also neglected. 

The general expression for the energy-momentum tensor that includes shear viscosity effects is 


T^ v — eu u u v — PA^ V + tt^ 


( 13 ) 
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where i is the shear stress tensor, u M is the fluid 4-velocity (u^u^ = 1), and the spatial projection operator is 
= 9iiv — u^Uy. We use the Landau definition for the local rest frame, = eidh The shear stress tensor 7r Mzy 

obeys a simplified Israel-Stewart relaxation equation 


*iuya(3 


Dn aP + -7 t^9 


+ Kuv = 277C7 


HV-) 


(14) 


W/3 




3 A/i^A, 


, the shear tensor 


O fiv — 


^a/ 3 V~"iP 4 the comoving derivative D = and tv = br]/(sT) is the shear viscosity relaxation time coefficient. 


where we have defined the tensor projector A Mm/ 3 = \ [A Ma 
A iivoLpS 01 " 1 ^, the comoving derivative D = and t. 

In this paper we take a const ant rj/ s that is adjusted in each case to obtain a reasonable description of LHC data 
(see the discussion in Section |iVB| ). This simplification overlooks the fact that r]/s may strongly depend on the 
temperature when T ~ 100 — 400 MeV [72H75] . which has been shown to affect elliptic flow m • We leave an analysis 
of the effects of a temperature dependent r]/s (and (/s) to a future publication. 

The v-USPhydro code accurate!}^ solves the energy-momentum conservation equations © (in their boost invari¬ 
ant form) with T^ v defined in © together with the relaxation equation ( [14] ) for i using the SPH algorithm. The 
number of SPH particles in the simulations carried out here ranged from 40000 - 50000 depending on the centrality 
class. The SPH h parameter for ideal hydrodynamics was set to be h = 0.1 fm while we used h = 0.3 fm in viscous 
simulations, as done in m- We used the lattice-based equation of state EOS S95n-vl from m and an isothermal 
Cooper-Frye [79? freezeout. The initial time for all the hydrodynamic simulations is To = 0.6 fm. Particle decays are 
taken into account using an adapted version of the AZHYDRO code [80! with hadronic resonances with masses up to 
1.7 GeV. Other technical details about our Lagrangian hydrodynamic solver can be found in m ■ 


A. Effects of the smoothing procedure on Knudsen and inverse Reynolds numbers 

The behavior displayed by the eccentricities indicates that the final azimuthal anisotropies, computed within viscous 
hydrodynamics, may be rather insensitive to variations in A at least in the mesoscopic regime. Before we present 
our results for the flow anisotropies, below we first show how the smoothing procedure we implemented change the 
Knudsen and inverse Reynolds numbers on an event by event basis. 

We compute the Knudsen number K u q = T n 6 [22] . which is the product of the relevant microscopic scale, tv, and 
a scalar that measures the spatial gradients of the flow, i.e, the expansion rate 9. We took the minimal value of tv in 
our viscous simulations, ~ 0.3 fm, as as lower bound on the value of the smoothing parameter A in viscous simulations 
to ensure that we are not evolving hydro dynamically truly microscopic sub-nucleon scales that are beyond the reach 
of Israel-Stewart hydro dynamic 

In Fig. [7] we show K nd = r n 6 within the radius of r < 5.5 fm for the same MCKLN event in Fig. [3] that belongs to 
the 0 — 5% centrality class of ^fs = 2.76 TeV Pb+Pb collisions at the LHC. On the left we set the smoothing scale 
of the initial conditions to be A = 0.3 fm (with r]/s = 0.11) and on the right A = 1 fm (where r]/s = 0.1125). The 
top panels were computed at r = 1 fm (i.e., 0.4 fm after the initial time) while the panels at the bottom show K u q 
at r = 4 fm. At early times (top panels), most of the matter is within the 5.5 fm radius we plotted and most of the 
regions where K u q > 1 form the edge of the system. However, note that even if one changes A by nearly an order of 
magnitude, the Knudsen number in the inner region is never small though most of details of the spatial structure are 
washed out when going from A = 0.3 to 1 fm. In fact, for A = 1 fm a disk of 2 fm radius centered at the origin would 
correspond to a nearly uniform K u q ^ 0.75 while when A = 0.3 fm the same 2 fm radius disk would still contain some 
subregions where K u q = 1. At later times (bottom panels), the system has had enough time to expand and the region 
where K u q > 1 has been pushed outside of the plotted disk, though the inner region still shows moderately large 
values around 0.5. These results show that the local Knudsen number decreases when one increases the smoothing 
scale of the initial condition (as expected), though the (lowest order) spatial eccentricities shown in Fig. [4] can remain 
nearly unchanged. 

To illustrate the issue with the applicability of hydrodynamics in event by event simulations containing initial state 
fluctuations of very short wavelength, we show in Fig. [8] the same event as before but now we take A = 0.1 fm. Even 
with the very small r]/s = 0.03 used in this calculation, one can see that the Knudsen number remains large and 
highly inhomogeneous throughout the hydrodynamic evolution. A comparison between Figs. [8] and [7] show that using 
the minimal value of r n ~ 0.3 fm in our simulations as a lower bound on A gives much more sensible Knudsen number 


13 We note that v-USPhydro was shown E] to reproduce both the analytical and semi-analytical radially expanding solutions of Israel- 
Stewart hydrodynamics first developed in E3- 

14 If one were to use the shear relaxation time computed at strong coupling in IZD which is roughly ~ 0.2/T, this lower bound would be 
approximately a factor of two smaller. This would essentially bring down the local K n g in Fig. □ by a factor of two. Therefore, one can 
see that the underlying assumptions about the microscopic nature of the fluid, i.e., weak (based on kinetic theory) or strong coupling, 
can shift up or down the ballpark estimate for the microscopic scale relevant for hydrodynamics. 
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FIG. 7. (Color online) Local Knudsen number, K u q = t n 0 , in the radius of r < 5.5 fm for the same MCKLN event in the 
0 — 5% centrality class of y/s — 2.76 TeV Pb+Pb collisions at the LHC displayed in Fig. [5] We use r n = brj/(sT) and initial 
time to = 0.6 fm. On the left we set A = 0.3 fm (with 77 /s = 0.11) and on the right A = 1 fm (with 77 /s = 0.1125). The top 
panels were computed at r = 1 fm while the panels at the bottom show K n e at r = 4 fm. 


profiles where the applicability of hydrodynamics can be at least suggested - the same cannot be said about the case 
where A = 0.1 fm. 

In Fig. [ 9 ] we show our results for the inverse Reynolds number, Re -1 = t^/P, for the same event. Once 

more, A = 0.3 fm on the left while on the right A = 1 fm. The top panels correspond to r = 1 fm and the bottom 
panels show Re -1 at r = 4 fm. Given that 7r M ^(ro,r) = 0, one can see that Re -1 is small in the 5.5 fm radius disk 
we plotted for both early and late times and different values of A. It would be interesting to investigate the case 
where 7 is not set to zero at To to see how its initial inhomogeneities and further time evolution are affected by the 
smoothing procedure considered in this paper. 


IV. RESULTS FOR THE ANISOTROPIC FLOW COEFFICIENTS IN NUCLEUS-NUCLEUS 

COLLISIONS 


We have seen in the previous sections that variations of the smoothing scale A by nearly an order of magnitude 
within the mesoscopic regime (see Fig. [I]) do not change significantly the initial spatial eccentricities in nucleus-nucleus 
collisions, though the local Knudsen number does decrease with increasing A. In this section we study how smoothing 
out the initial conditions affect the anisotropic flow coefficients in Pb+Pb collisions. We use the event plane method 
j8Tj to compute the anisotropic flow coefficients v n (see Appendix C of [40. for details) of the 150 individual events in 
each centrality class. 
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FIG. 8. (Color online) Local Knudsen number, K n o — tv 0 , in the radius of r < 5.5 fm for the same MCKLN event in the 
0 — 5% centrality class of y/s = 2.76 TeV Pb+Pb collisions at the LHC displayed in Fig. [5] We use again = 5rj/(sT) and 
initial time to = 0.6 fm though now A = 0.1 fm and rj/s = 0.03. On the left show K n e at r = 1 fm while on the right r — 4 fm. 


A. Flow coefficients within ideal hydrodynamics 


Here we study how changes in A affect the flow harmonics computed in the inviscid limit, i.e., r]/s = 0. After 
performing the smoothing procedure on the initial energy density, the particle spectrum (before integration over the 
azimuthal angle) acquires a dependence on the smoothing scale A, which is later passed over to the final u n ’s obtained 
after freezeout+decay. We have set the freezeout temperature to be Tpo = 130 MeV in the ideal hydrodynamic 
calculations. 

In Fig. 10 we show how the spectrum dN/(‘hrpTdpTdy) of all charged particles for the 0 — 5% centrality class, 
computed using ideal hydrodynamics, is affected when A varies from 0.1 to 1 fm. CMS data [82], 83]] is included in 
this plot to give an idea about the magnitude of the effect. While at sufficiently low pr < 0.5 GeV the spectrum can 
be made robust with respect to variations in A (by appropriately choosing the overall constant in the initial energy 
density), as one increases A the high p T part is depleted and the spectrum gets softer. This is expected since the short 
wavelength fluctuations present in event by event simulations that generally contribute to the high pr part of the 
spectrum [61] 63; are systematically removed when going from A = 0.1 to 1 fm, which in turn decreases the average 
transverse momentum. _ 

In Fig.llllwe show how the differential flow harmonics of all charged particles v u (pt) = yj(v^ipr)} (computed within 
ideal hydrodynamics and averaged over 150 events in each centrality class) change as one increases the smoothing 
parameter A from 0.1 to 1 fm. The left panels correspond to the 0 — 5% centrality class while the right panels show 
the 20 — 30% centrality results for y/s = 2.76 TeV Pb+Pb collisions at the LHC. Once more, the data points here 
correspond to CMS data [82j |83], which are only included at this point to give an idea of the magnitude of the 
generated anisotropic flow. There are a few salient features encoded in these plots. First, in accordance with the 
results for the spectrum discussed before, one sees that the anisotropic flow coefficients are enhanced when short 
wavelength fluctuations are smoothed outpq In our case, the addition of short wavelength fluctuations bring the u n ’s 
down because their contribution to the overall anisotropy is not sufficient to overcome the increase in the particle 
spectrum in Fig. [l0]for small values of A. Note also that ideal hydrodynamics with A = 0.1 fm is already very close 
to the data (vs is actually below the data) especially for 0 — 5% centrality. Another interesting feature displayed in 
Fig. 11 is that v u (pt) changes less when A goes from 0.3 to 1 fm than it does when A varies between 0.1 and 0.3 fm. 
Moreover, we found that odd harmonics converge faster than their even counterparts. 

In Fig. [12] we show the corresponding results for the px-integrated flow harmonics of all charged particles, v n = 
yj(v^) (computed within ideal hydrodynamics and averaged over 150 events in each centrality class), and their 
dependence on A. The left panel corresponds to the 0 — 5% centrality class while the right panel shows the 20 — 30% 
centrality results for y/s = 2.76 TeV/n Pb+Pb collisions at the LHC. For the most central collisions, the robustness of 
the initial eccentricities in Figs. [4] and [5] are transferred to the u n ’s, which are essentially flat with respect to variations 


15 Ref. j62t observed the same behavior in their calculations of elliptic flow. 
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FIG. 9. (Color online) Local inverse Reynolds number, Re -1 m in the radius of r < 5.5 fm for the same MCKLN 

event in the 0 — 5% centrality class of y/s — 2.76 TeV Pb+Pb collisions at the LHC displayed in Fig. [ 3 ] We use tv = 5tj/(sT) 
and initial time To = 0.6 fm at which 7 +^ = 0. On the left we set the smoothing scale to be A = 0.3 fm (with 77 /s = 0.11) and 
on the right A = 1 fm (with 77 /s = 0.1125). The top panels were computed at t = 1 fm while the panels at the bottom show 
Re~ 1 at t = 4 fm. 


in A from 0.1 to 1 fm. On the other hand, for more peripheral collisions (left panel) we find an enhancement in v 2 
for larger values of A (this effect is considerably reduced when shear viscosity is taken into account, see IVB). This 
indicates that, at least in the case of MCKLN initial conditions, dv n (X)/d\ > 0 and the system’s response to the initial 
spatial eccentricities is somewhat improved when very short wavelength fluctuations are removed and only scales in 
the mesoscopic regime in Fig. [l]are taken into account. It would be interesting to check if this statement holds for 
other initial state models as well or if it depends on the specific origin of the initial state fluctuations. 


B. Flow coefficients within viscous hydrodynamics and comparison to LHC data 

In this section we show our viscous hydrodynamic results for A = 0.3 fm and A = 1 fm. In order to better 
describe CMS data, we used 77/s = 0.11 for A = 0.3 fm while for A = 1 fm we used 77/s = 0.1125. In these viscous 
simulations, the freezeout temperature is Tpo = 120 MeV (unless stated otherwise) and the shear viscous contribution 
to the distribution function used in the Cooper-Frye freezeout is Sf = 7r /iu p^'p l/ /(2sT 3 ) for all hadronic species. The 
spectrum of all charged particles computed in the 0 — 5% centrality class for yfs = 2.76 TeV Pb+Pb collisions at 
the LHC can be seen in Fig. [T3| where we also included the A = 0.1 fm ideal hydrodynamic curve from Fig. [To| for 
comparison. As before, the data points correspond to CMS data for all charged particles [82, 83] • One can see that 
the nonzero shear viscosity did not change the qualitative behavior seen already in the ideal case in Fig. [lOj the low 
Pt yield is barely affected when one increases A and (small) differences can only be seen above p T = 1 GeV. 
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Pt [GeV] 


FIG. 10. (Color online) dN/(2jvpTdpTdy) of all charged particles in the 0 — 5% centrality class of /s = 2.76 TeV Pb+Pb 
collisions at the LHC computed using ideal hydrodynamics (averaged over 150 events) for A = 0.1 — 1.0 fm. Data points 
correspond to CMS data [52] [83 . 


The results for the differential flow harmonics V 2 (jpr) and v%( jpt) f° r the 20 — 30% centrality class, computed now 
within viscous hydrodynamics, can be found in Fig. [l4j One can see that once shear viscosity is included there is 
almost no difference between the curves computed using A = 0.3 fm from those where A = 1 fm. Both values of A give 
a good description of the CMS data even though the initial energy densities and Knudsen numbers corresponding to 
these curves (see Figs. [3] and [7]) reflect very different levels of spatial resolution. 

In Fig. [l5] we show the corresponding results for the pr integrated v/s of all charged particles in the 20 — 30% 
centrality class. The grey bars correspond to CMS data [82] [83] with the bottom bars denoting the 20 — 25% centrality 
while the top bars stand for the 25 — 30% centrality class. One can see that the integrated v/s computed within event- 
by-event viscous hydrodynamics change very little when A is in the mesoscopic regime. The small error bands attached 
to the red and black points show how an increase in the freezeout temperature from 120 MeV to 130 MeV affects the 
u n ’s (i.e., they decrease if Tpo increases). Also, the viscous hydrodynamic calculations for the flow anisotropies are 
in the ballpark of the CMS data except for vs though the ideal hydrodynamic calculation computed using A = 0.1 fm 
(the blue points included in the plot for comparison) is below the data for all n. 

The results of this section give support to the idea that viscous hydrodynamic calculations of anisotropic flow 
coefficients are remarkably robust with respect to the spatial resolution scale present in the initial conditions. This 
shows that the short wavelength fluctuations present in MCKLN initial conditions, which lead to a large Knudsen 
number in event by event simulations, can be consistently smoothed out without changing significantly the azimuthal 
anisotropies responsible for generating the flow coefficients v n . It would be interesting to check how the v n distri¬ 
butions computed within hydrodynamics [84] |85] are affected by variations in A. Such a study requires performing 
hydrodynamic calculations in a number of events an order of magnitude larger than the one presented here and is 
beyond the scope of this paper. 


V. SMOOTHING OUT PROTON-NUCLEUS COLLISIONS 


We saw in the previous sections that the robustness of the initial eccentricities with respect to variations of the 
scale of spatial fluctuations in the mesoscopic regime in nucleus-nucleus collisions (across centrality classes) is carried 
over to the final anisotropic flow coefficients. This is the expected outcome based on the discussion in [84] • However, 
with the advent of the high multiplicity p+Pb data at the LHC [53] and the subsequent observation of significant 
collective multi-particle behavior in these collisions |86| , it is interesting to investigate how sensitive the eccentricities 
and the final flow harmonics in p+Pb collisions are to spatial fluctuations at small scales El- 

In this paper we use the Trento code [52] to generate top 1% high multiplicity /s = 5.02 TeV p+Pb events and 
check how the eccentricities computed using this model vary with the smoothing parameter. We fix the parameters 
of the model to be p = 0.3 and k = 1.3, which were shown to provide a good description of the number of charged 
particles per unit rapidity found in different collision systems and centrality classes [52] . We show on the left panel 
in Fig. 16 our results for the A dependence of £ 2 , n corresponding to the top 1% high multiplicity /s = 5.02 TeV 
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FIG. 11. (Color online) Convergence of the differential flow harmonics of all charged particles v n (pr) = \/Jv?JjPt)) (computed 
within ideal hydrodynamics and averaged over 150 events in each centrality class) as one increases the smoothing parameter A 
from 0.1 to 1 fm. The left panels correspond to the 0 — 5% centrality class while the right panels show the 20 — 30% centrality 
results for = 2.76 TeV Pb+Pb collisions at the LHC. The data points correspond to CMS data [82] 155]. 


p+Pb collisions computed using Trento. On the right panel we show the respective MCKLN results for the 65 — 70% 
centrality class of yfs = 2.75 TeV Pb+Pb collisions, which have comparable multiplicities [53] to the p+Pb events. 
One can see that the p+Pb system is much more sensitive to mesoscopic scale fluctuations than peripheral Pb+Pb, 
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FIG. 12 . (Color online) Dependence of the integrated flow harmonics v n = yj (v%) (computed within ideal hydrodynamics and 
averaged over 150 events in each centrality class) on the smoothing scale A. The left panel corresponds to the 0 — 5% centrality 
class while the right panel shows the 20 — 30% centrality results for y/s = 2.76 TeV Pb+Pb collisions at the LHC. 




FIG. 13. (Color online) dN/ {2'Kprdprdy) of all charged particles in the 0 — 5% centrality class of y/s = 2.76 TeV Pb+Pb collisions 
at the LHC computed using shear viscous hydrodynamics (averaged over 150 events) with A = 0.3 fm (and rj/s = 0.11) and 
A = 1 fm (where 77 /s = 0.1125). The previous ideal hydrodynamic curve computed with A = 0.1 fm is also included for 
comparison. Data points correspond to CMS data [82] [83]. 


even at the same multiplicity. This is to be expected given the smaller size of p+Pb and that in this system the 
separation of scales in Fig. [l]does not apply. Another interesting feature observed in these plots is that the overall 
magnitude of the eccentricities in p+Pb is a factor of two smaller than those found in peripheral Pb+Pb at the same 
multiplicity (we note that the p+Pb eccentricities we found using the Trento code are compatible with those found in 
IP-Glasma [45]). This puts into perspective the experimental finding [53] that the two particle cumulant v% is nearly 
the same in these systems. 

The fact that p+Pb is much more sensitive to the smoothing parameter than peripheral Pb+Pb can be more 
clearly seen in Fig. [17] where we show the relative variation of the eccentricities with A. While in peripheral Pb+Pb 
the variation of the eccentricities remain below 10% for A = 0.1 — 1 fm, the same cannot be said about p+Pb where 
10% variations are already found in the interval A = 0.1 — 0.3 fm (note, however, that the variations in the lowest 
order eccentricities £2,2 and £2,3 remain ^ 10% up to A = 0.5 fm). We note here that if one uses the Trento setup 
specific to LHC found in [52] of p = 0 and k = 1.4 that the initial eccentricities are even smaller and even more 
sensitive to the smoothing scale and even see large deviations for £2,2 and £2,3 above A = 0.3. This shows that p+Pb 
collisions should be very sensitive to the underlying physics below the confinement scale, which may justify the search 
for alternative non-hydrodynamical explanations for the p+Pb flow harmonics data [88H94] . 
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FIG. 14. (Color online) Differential flow harmonics V2(pt) (top, solid curves) and V3(pr) (bottom, solid curves) of all charged 
particles in the 20 — 30% centrality class for y/s = 2.76 TeV Pb+Pb collisions at the LHC. The left panels were computed using 
A = 0.3 fm while the right panels used A = 1 fm in the calculations. The data points correspond to CMS data [82j [83] and the 
solid lines are the shear viscous results for Tfo — 120 MeV. Ideal hydrodynamic results for the flow coefficients from Fig. mi 
are shown in dashed lines for comparison. 


VI. CONCLUSIONS AND OUTLOOK 

We argued in this paper that QCD phenomena in ultrarelativistic collisions of heavy ions may be separated into 
three distinct regimes, illustrated in Fig. [l] The microscopic sub-nucleon scale regime involves phenomena that are 
sensitive to length scales of the order of 1/Q S while one may take 1/Aqcd as a mesoscopic scale above which one 
encounters a macroscopic nuclear scale regime given by length scales of the order of the radius of a large nucleus. 
Though there is currently no single effective theory that is able to describe all the different phenomena that occur 
at these different length scales, the current success of viscous hydrodynamics in understanding several features of 
heavy ion data implies that it is should be a good starting point to describe the bulk dynamics in the mesoscopic 
regime since that is actually probed in event by event calculations. In this case, it is important to remark that current 
viscous hydrodynamic implementations based on Israel-Stewart hydrodynamics do not have the correct degrees of 
freedom to describe scales shorter than the respective relaxation time (in our case, r^), which provides a lower bound 
on the length scales at which these hydrodynamic models are sensible. The presence of these different scales in the 
description of ultrarelativistic heavy ion collisions calls for a systematic study of the effects of initial state fluctuations, 
separated according to their wavelength [46], on hydro dynamically generated observables such as the anisotropic flow 
coefficients. 

In this paper we used a Lagrangian based procedure to smooth out the initial conditions employed in hydrodynamics 
simulations of ultrarelativistic heavy ion collisions that allowed us to investigate, in a systematic manner, how the final 
azimuthal anisotropies u n ’s (computed within event by event hydrodynamics) depend on the wavelength of the initial 
energy density fluctuations. The smoothing procedure is characterized by a smearing function with a finite support 
given by the smoothing scale A that defines the maximal size of the initial energy density gradients. The method 
can be applied in a variety of initial condition models ranging from MCKLN/MCGlauber to IP-Glasma. Here we 
used y/s = 2.76 TeV Pb+Pb MCKLN initial conditions to show that the initial spatial eccentricities that govern the 
final state flow harmonics are very robust with respect to variations in the underlying scale of initial energy density 
fluctuations (or smoothing parameter) A, as long as these fluctuations are in the so-called mesoscopic regime (see Fig. 
[l]) where 0.1 < A < 1 fm. 

Given that the local Knudsen number is K n ~ 1/A, the robustness of the initial eccentricities with respect to changes 
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FIG. 15. (Color online) Dependence of the integrated flow harmonics v n of all charged particles (computed within viscous 
hydrodynamics and averaged over 150 events) on the smoothing scale A. The black points were computed using A = 0.3 fm 
(with rj/s = 0.11) while for the red points A = 1 fm (and rj/s = 0.1125). The small error band attached to these points denote 
the small dependence of the results when one increases the freezeout temperature from 120 MeV to 130 MeV. The grey bars 
correspond to CMS data [82], 83] for y/s = 2.76 TeV Pb+Pb collisions at the LHC with the bottom bars denoting the 20 — 25% 
centrality while the top bars stand for the 25 — 30% centrality class (note that for V 3 , V 4 , and V 5 the bars overlap). We have 
included the ideal hydrodynamic result for A = 0.1 fm from Fig. [12] for a comparison. 
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FIG. 16. (Color online) Dependence of the initial eccentricities £ 2 , n = y (e|, n ) on fh e smoothing parameter, A, for the top 1% 

high multiplicity y/s — 5.02 TeV p+Pb collisions (left panel) compared to the same quantity in 65 — 70% centrality class of 
y/s = 2.76 TeV Pb+Pb events (right panel). These systems have comparable multiplicities [53]. The average in each case was 
performed over 150 events. 


in the fluctuation scale in the mesoscopic regime was carried over to the r+’s computed within viscous hydrodynamics. 
This implies that events with large local K n and events where K n is near the hydrodynamic regime lead to nearly 
indistinguishable flow harmonics. In fact, at least in the case of Pb+Pb MCKLN initial conditions (and most certainly 
also for MCGlauber), the anisotropic flow coefficients computed within event by event viscous hydrodynamics are only 
sensitive to long wavelength scales of order 1/A qcd ^ 1 fm and are, thus, very robust with respect to variations in 
the initial local Knudsen number. 

Since the structures at very small scales present in MCKLN initial conditions are still invariably rooted in the 
MC Glauber-like fluctuations of the positions of the (independent) nucleons in the nucleus, one may say that 1/A qcd 
should indeed be the upper limit of the mesoscopic regime in our case. Therefore, it is important to investigate if the 
conclusions found here regarding the robustness of the flow harmonics in viscous hydrodynamics in nucleus-nucleus 
collisions still hold when the short wavelength fluctuations in the initial energy density stem from nontrivial color 
fields defined in the microscopic sub-nucleon scale regime ~ 1/Q S shown in Fig. [l] One should keep in mind, however, 
that taking into account fluctuations of very short wavelength scales considerably worsens the local Knudsen number 
(see Fig. [8|. Moreover, since the relaxation time provides a lower limit on the scales that can be described within 
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FIG. 17. (Color online) Relative variation of initial eccentricities A £ 2 ,n = 100 [£2,n(Ao) — £2,n(A)] /£2,n(Ao) on the smoothing 
parameter, A, for the top 1% high multiplicity y / s = 5.02 TeV p+Pb collisions (left panel) compared to the same quantity in 
65 — 70% centrality class of yfs = 2.76 TeV Pb+Pb events (right panel). The average in each case was performed over 150 
events. 


Israel-Stewart hydrodynamics, in order to assess very short wavelength fluctuations characterized by A < 0.1 fm, r^ 
would have to be significantly small. 


Furthermore, effective theory arguments may be used to argue that truly microscopic sub-nucleon ~ 1/Q S physics 
should not influence the bulk dynamics of matter at much larger length scales ^ 1/A qcd, which seem to be the 
relevant ones for anisotropic flow coefficients in nucleus-nucleus collisions where the separation of scales shown in Fig. 
[l] approximately holds. However, it is known that in many-body systems sometimes the existence of a large separation 
of scales is not sufficient to rule out microscopic effectand, thus, it would definitely be interesting to check how the 


anisotropic flow generated with other type of initial conditions is affected by the smoothing procedure defined in |II A 
Also, it would be interesting to understand how fluctuations in the sub-nucleon regime may affect the approximate 
scaling properties found in heavy ion observables discussed in da nz] and also how they can influence a principal 
component analysis [981199] . Furthermore, we intend to investigate the case of nucleus-nucleus collisions with reduced 
center of mass collision energies y/s, such as in the low energy beam energy scan program at RHIC. 


We also investigated how the eccentricities of the top 1% high multiplicity p+Pb collisions are affected by our 
smoothing procedure. In this case, the scenario depicted in Fig. [l] does not hold since the separation between the 
macroscopic and the mesoscopic nuclear regimes becomes artificial and much less motivated than in the case of nucleus- 
nucleus collisions. We found that the eccentricities in p+Pb are very sensitive to sub-nucleon scale fluctuations, which 
should be contrasted with the robustness found in peripheral Pb+Pb collisions with the same multiplicity. 


Also, regarding the hydrodynamic evolution, it is important to take into account the temperature dependence of 
transport coefficients such as r]/s (and (/s ) since that will affect our results for the Knudsen number [22]. The 
inclusion of other second order hydrodynamic coefficients beyond in the hydrodynamic evolution may also play a 
role in determining iT n , especially in small systems. Moreover, it would be interesting to extend the study performed 
here to the case of 3+1 viscous hydrodynamic simulations including the possibility of different smoothing scales in 
the transverse and longitudinal directions. 

Since the underlying bulk dynamics of the expanding medium is relevant for energy loss calculations [9], it would 
be interesting to check to what extent jet quenching-related observables are sensitive to short wavelength fluctuations. 
An initial study performed in m showed little sensitivity to initial state fluctuations, though it should be kept in 
mind that a full event-by-event hydrodynamic evolution of the initial spatial inhomogeneities was not included in 
their analysis. Additionally, one may also study how recent calculations of the medium response to the energy and 
momentum lost by jets [10114103] are affected by different levels of spatial resolution parametrized in terms of the 
smoothing scale A. 


16 The stability of matter (see [95]) relies on the Pauli principle, which is not only essential in the description of matter in subatomic scales 
but is also ultimately fundamental to ensure the stability of macroscopic bodies. 












20 


ACKNOWLEDGMENTS 


We thank T. Kodama for general discussions about the coarse-graining process in heavy ion collisions, A. Du- 
mitru and G. S. Denicol for enlightening discussions about the effects of sub-nucleon initial state fluctuations in 
hydrodynamics, G. Torrieri for suggesting to check the role of fluctuations in different collision systems and energies, 
M. Luzum for discussions on the hydrodynamic approach to proton-nucleus collisions, and S. Moreland and J. Bern- 
hard for assistance with using the Trento initial condition code. JNH and MG acknowledge support from the US-DOE 
Nuclear Science Grant No. DE-FG02-93ER40764. JN thanks the Physics Department at Columbia University for 
its hospitality and Fundacao de Amparo a Pesquisa do Estado de Sao Paulo (FAPESP) and Conselho Nacional de 
Desenvolvimento Cientffico e Tecnologico (CNPq) for support. 


[1] BRAHMS collaboration, I. Arsene et al., Nucl. Phys. A 757, 1 (2005), arXiv:nucl-ex/0410020 . 

[2] PHENIX collaboration, K. Adcox et al., Nucl. Phys. A 757, 184 (2005), arXiv:nucl-ex/0410003 . 

[3] PHOBOS collaboration, B. B. Back et al., Nucl. Phys. A 757, 28 (2005), arXiv:nucl-ex/0410022 . 

[4] STAR collaboration, J. Adams et al., Nucl. Phys. A 757, 102 (2005), arXiv:nucl-ex/0501009 . 

[5] M. Gyulassy and L. McLerran, Nucl. Phys. A 750, 30 (2005), arXiv:nucl-th/0405013 . 

[6] J. D. Bjorken, FERMILAB-PUB-82-059-THY, FERMILAB-PUB-82-059-T. 

[7] M. Gyulassy and M. Plumer, Phys. Lett. B 243, 432 (1990). 

[8] X. N. Wang and M. Gyulassy, Phys. Rev. Lett. 68, 1480 (1992). 

[9] K. M. Burke et al. [JET Collaboration], Phys. Rev. C 90, no. 1, 014909 (2014) arXiv: 1312.5003 [nucl-th]]. 

[10] U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013), arXiv:1301.2826 [nucl-th]]. 

[11] P. Mota, T. Kodama, R. Derradi de Souza and J. Takahashi, Eur. Phys. J. A 48, 165 (2012) arXiv: 1210.3129 [hep-ph]]. 

[12] R. D. de Souza, T. Koide and T. Kodama, arXiv: 1506.03863 [nucl-th]. 

[13] S. Weinberg, “Cosmology”, Oxford University Press, 2008. 

[14] F. Gelis, E. Iancu, J. Jalilian-Marian and R. Venugopalan, Ann. Rev. Nucl. Part. Sci. 60, 463 (2010) arXiv:1002.0333 
[hep-ph]]. 

[15] W. Israel and J. M. Stewart, Annals Phys. 118, 341 (1979). 

[16] G. S. Denicol, T. Koide and D. H. Rischke, Phys. Rev. Lett. 105, 162501 (2010) arXiv: 1004.5013 [nucl-th]]. 

[17] G. S. Denicol, J. Noronha, H. Niemi and D. H. Rischke, Phys. Rev. D 83, 074019 (2011) arXiv:1102.4780 [hep-th]]. 

[18] G. S. Denicol, H. Niemi, E. Molnar and D. H. Rischke, Phys. Rev. D 85, 114047 (2012) [Phys. Rev. D 91, no. 3, 039902 
(2015)] arXiv:1202.4551 [nucl-th]]. 

[19] J. Noronha and G. S. Denicol, arXiv:1104.2415 [hep-th]. 

[20] M. P. Heller, R. A. Janik, M. Spali?ski and P. Witaszczyk, Phys. Rev. Lett. 113, no. 26, 261601 (2014) arXiv: 1409.5087 
[hep-th]]. 

[21] M. Nopoush, M. Strickland, R. Ryblewski, D. Bazow, U. Heinz and M. Martinez, arXiv: 1506.05278 [nucl-th]. 

[22] H. Niemi and G. S. Denicol, arXiv: 1404.7327 [nucl-th]. 

[23] J. Casalderrey-Solana, E. V. Shuryak and D. Teaney, J. Phys. Conf. Ser. 27, 22 (2005) [Nucl. Phys. A 774, 577 (2006)] 
hcp-ph/0411315 . 

[24] A. Yarom, Phys. Rev. D 75, 105023 (2007) hep-th/0703095 . 

[25] P. M. Chesler and L. G. Yaffe, Phys. Rev. Lett. 99, 152001 (2007) arXiv:0706.0368 [hep-th]]. 

[26] S. S. Gubser, S. S. Pufu and A. Yarom, Phys. Rev. Lett. 100, 012301 (2008) arXiv:0706.4307 [hep-th]]. 

[27] J. Noronha, G. Torrieri and M. Gyulassy, Phys. Rev. C 78, 024903 (2008) arXiv:0712.1053 [hep-ph]]. 

[28] J. Noronha, M. Gyulassy and G. Torrieri, Phys. Rev. Lett. 102, 102301 (2009) arXiv:0807.1038 [hep-ph]]. 

[29] R. B. Neufeld, B. Muller and J. Ruppert, Phys. Rev. C 78, 041901 (2008) arXiv:0802.2254 [hep-ph]]. 

[30] R. B. Neufeld, Phys. Rev. D 78, 085015 (2008) arXiv:0805.0385 [hep-ph]]. 

[31] B. Betz, M. Gyulassy, J. Noronha and G. Torrieri, Phys. Lett. B 675, 340 (2009) arXiv:0807.4526 [hep-ph]]. 

[32] B. Betz, J. Noronha, G. Torrieri, M. Gyulassy, I. Mishustin and D. H. Rischke, Phys. Rev. C 79, 034902 (2009) 
arXiv:0812.4401 [nucl-th]]. 

[33] M. L. Miller, K. Reygers, S. J. Sanders and P. Steinberg, Ann. Rev. Nucl. Part. Sci. 57, 205 (2007) nucl-ex/0701025 . 

[34] B. Alver, M. Baker, C. Loizides and P. Steinberg, arXiv:0805.4411 [nucl-ex]. 

[35] C. Loizides, J. Nagle and P. Steinberg, arXiv: 1408.2549 [nucl-ex]. 

[36] H.-J. Drescher and Y. Nara, Phys. Rev. C 75, 034905 (2007) nucl-th/0611017 . 

[37] B. Schenke, P. Tribedy and R. Venugopalan, Phys. Rev. Lett. 108, 252301 (2012) arXiv: 1202.6646 [nucl-th]]. 

[38] A. Dumitru and Y. Nara, Phys. Rev. C 85, 034907 (2012) arXiv: 1201.6382 [nucl-th]]. 

[39] J. Noronha-Hostler, J. Noronha, G. S. Denicol, R. P. G. Andrade, F. Grassi and C. Greiner, Phys. Rev. C 89, no. 5, 
054904 (2014) arXiv: 1302.7038 [nucl-th]]. 

[40] J. Noronha-Hostler, G. S. Denicol, J. Noronha, R. P. G. Andrade and F. Grassi, Phys. Rev. C 88, 044916 (2013) 
arXiv: 1305.1981 [nucl-th]]. 

[41] J. Noronha-Hostler, J. Noronha and F. Grassi, Phys. Rev. C 90, no. 3, 034907 (2014) arXiv: 1406.3333 [nucl-th]]. 





21 


[42] H. Petersen, C. Coleman-Smith, S. A. Bass and R. Wolpert, J. Phys. G 38, 045102 (2011) arXiv: 1012.4629 [nucl-th]]. 

[43] C. E. Coleman-Smith, H. Petersen and R. L. Wolpert, J. Phys. G 40, 095103 (2013) arXiv: 1204.5774 [hep-ph]]. 

[44] M. Rihan Haque, V. Roy and A. K. Chaudhuri, Phys. Rev. C 86, 037901 (2012) arXiv: 1204.2986 [nucl-ex]]. 

[45] A. Bzdak, B. Schenke, P. Tribedy and R. Venugopalan, Phys. Rev. C 87, no. 6, 064906 (2013) arXiv: 1304.3403 [nucl-th]]. 

[46] S. Floerchinger and U. A. Wiedemann, Phys. Lett. B 728, 407 (2014) arXiv: 1307.3453 [hep-ph]]. 

[47] S. Floerchinger and U. A. Wiedemann, Phys. Rev. C 88, 044906 (2013) arXiv:1307.7611 [hep-ph]]. 

[48] S. Floerchinger and U. A. Wiedemann, Phys. Rev. C 89, no. 3, 034914 (2014) arXiv:1311.7613 [hep-ph]]. 

[49] E. Retinskaya, M. Luzum and J. Y. Ollitrault, Phys. Rev. C 89, no. 1, 014902 (2014) arXiv:1311.5339 [nucl-th]]. 

[50] T. Renk and H. Niemi, Phys. Rev. C 89, no. 6, 064907 (2014) arXiv: 1401.2069 [nucl-th]]. 

[51] V. P. Konchakovski, W. Cassing and V. D. Toneev, J. Phys. G 42, no. 5, 055106 (2015) arXiv: 1411.5534 [nucl-th]]. 

[52] J. S. Moreland, J. E. Bernhard and S. A. Bass, Phys. Rev. C 92, no. 1, 011901 (2015) arXiv:1412.4708 [nucl-th]]. 

[53] S. Chatrchyan et al. [CMS Collaboration], Phys. Lett. B 724, 213 (2013) arXiv: 1305.0609 [nucl-ex]]. 

[54] Richard S. Hamilton, “Three-manifolds with positive Ricci curvature,” J. Differential Geom. 17, no. 2, 255-306 (1982). For 
an introduction see P. Topping, Lectures on the Ricci Flow , London Mathematical Society Lecture Note Series, Cambridge 
University Press, UK 2006. 

[55] L. B. Lucy, ApJ. 82, 1013 (1977); R. A. Gingold and J. J. Monaghan, MNRAS 181, 375 (1977); R. A. Gingold and 
J. J. Monaghan, Mon. Not. Roy. Astron. Soc. 181, 375 (1977). 

[56] J. J. Monaghan, Amur. Rev. Astron. Astrophys. 30, 543 (1992); E. Chow and J. J. Monaghan, J. Comput. Phys. 134, 296 
(1997); P. Laguna, W. A. Miller and W. H. Zurek, Phys. Rev. D 41, 451 (1990); R. J. Thacker, E. R. Tittley, F. R. Pearce, 
H. M. P. Couchman and P. A. Thomas, Mon. Not. Roy. Astron. Soc. 319, 619 (2000). 

[57] C. E. Aguiar, T. Kodama, T. Osada and Y. Hama, J. Phys. G 27, 75 (2001). 

[58] T. Osada, C. E. Aguiar, Y. Hama and T. Kodama, nucl-th/0102011 

[59] C. E. Aguiar, Y. Hama, T. Kodama and T. Osada, Nucl. Phys. A 698, 639 (2002). 

[60] O. Socolowski, Jr., F. Grassi, Y. Hama and T. Kodama, Phys. Rev. Lett. 93, 182301 (2004). 

[61] Y. Hama, T. Kodama and O. Socolowski, Jr., Braz. J. Phys. 35, 24 (2005). 

[62] R. Andrade, F. Grassi, Y. Hama, T. Kodama and O. Socolowski, Jr., Phys. Rev. Lett. 97, 202302 (2006). 

[63] R. P. G. Andrade, F. Grassi, Y. Hama, T. Kodama and W. L. Qian, Phys. Rev. Lett. 101, 112301 (2008). 

[64] J. Takahashi, B. M. Tavares, W. L. Qian, R. Andrade, F. Grassi, Y. Hama, T. Kodama and N. Xu, Phys. Rev. Lett. 103, 
242301 (2009). 

[65] F. G. Gardim, F. Grassi, M. Luzum and J. -Y. Ollitrault, Phys. Rev. Lett. 109, 202302 (2012). 

[66] R. S. Bhalerao, M. Luzum and J. Y. Ollitrault, Phys. Rev. C 84, 054901 (2011) arXiv: 1107.5485 [nucl-th]]. 

[67] F. G. Gardim, F. Grassi, M. Luzum and J. Y. Ollitrault, Phys. Rev. C 85, 024908 (2012) arXiv:1111.6538 [nucl-th]]. 

[68] F. G. Gardim, J. Noronha-Hostler, M. Luzum and F. Grassi, Phys. Rev. C 91, no. 3, 034902 (2015) arXiv: 1411.2574 

[nucl-th]]. 

[69] J. B. Rose, J. F. Paquet, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon and C. Gale, Nucl. Phys. A 931, 926 (2014) 
arXiv: 1408.0024 [nucl-th]]. 

[70] S. Ryu, J.-F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon and C. Gale, arXiv: 1502.01675 [nucl-th]. 

[71] S. I. Finazzo, R. Rougemont, H. Marrochio and J. Noronha, JHEP 1502, 051 (2015) arXiv:1412.2968 [hep-ph]]. 

[72] T. Hirano and M. Gyulassy, Nucl. Phys. A 769, 71 (2006) nucl-th/0506049 . 

[73] L. P. Csernai, J. I. Kapusta and L. D. McLerran, Phys. Rev. Lett. 97, 152303 (2006) nucl-th/0604032 . 

[74] J. Noronha-Hostler, J. Noronha and C. Greiner, Phys. Rev. Lett. 103, 172302 (2009) arXiv:0811.1571 [nucl-th]]. 

[75] J. Noronha-Hostler, J. Noronha and C. Greiner, Phys. Rev. C 86, 024913 (2012) arXiv:1206.5138 [nucl-th]]. 

[76] H. Niemi, G. S. Denicol, P. Huovinen, E. Molnar and D. H. Rischke, Phys. Rev. Lett. 106, 212302 (2011) arXiv: 1101.2442 
[nucl-th]]. 

[77] H. Marrochio, J. Noronha, G. S. Denicol, M. Luzum, S. Jeon and C. Gale, Phys. Rev. C 91, no. 1, 014903 (2015) 
arXiv: 1307.6130 [nucl-th]]. 

[78] P. Huovinen and P. Petreczky, Nucl. Phys. A 837, 26 (2010) arXiv:0912.2541 [hep-ph]]. 

[79] F. Cooper and G. Frye, Phys. Rev. D 10, 186 (1974). 

[80] P. F. Kolb, J. Sollfrank, and U. Heinz, Phys. Rev. C 62 (2000) 054909; P. F. Kolb and R. Rapp, Phys. Rev. C 67 (2003) 
044903; P. F. Kolb and U. Heinz, nucl-th/0305084 

[81] A. M. Poskanzer and S. A. Voloshin, Phys. Rev. C 58, 1671 (1 

[82] S. Chatrchyan et al. [CMS Collaboration], Phys. Rev. C 87, nc 

[83] S. Chatrchyan et al. [CMS Collaboration], Phys. Rev. C 89, nc 

[84] H. Niemi, G. S. Denicol, H. Holopainen and P. Huovinen, F 

[nucl-th]]. 

[85] H. Niemi, K. J. Eskola and R. Paatelainen, arXiv: 1505.02677 [hep-ph]. 

[86] V. Khachatryan et al. [CMS Collaboration], Phys. Rev. Lett. 115, no. 1, 012301 (2015) arXiv:1502.05382 [nucl-ex]]. 

[87] B. Schenke and R. Venugopalan, Phys. Rev. Lett. 113, 102301 (2014) arXiv: 1405.3605 [nucl-th]]. 

[88] K. Dusling and R. Venugopalan, Phys. Rev. D 87, no. 5, 054014 (2013) arXiv: 1211.3701 [hep-ph]]. 

J. Noronha and A. Dumitru, Phys. Rev. D 89, no. 9, 094008 (2014) arXiv: 1401.4467 [hep-ph]]. 

M. Gyulassy, P. Levai, I. Vitev and T. S. Biro, Phys. Rev. D 90, no. 5, 054025 (2014) arXiv: 1405.7825 [hep-ph]]. 

A. Dumitru and A. V. Giannini, Nucl. Phys. A 933, 212 (2014) arXiv: 1406.5781 [hep-ph]]. 

A. Dumitru, L. McLerran and V. Skokov, Phys. Lett. B 743, 134 (2015) arXiv: 1410.4844 [hep-ph]]. 

T. Lappi, Phys. Lett. B 744, 315 (2015) arXiv:1501.05505 [hep-ph]]. 


nucl-ex/9805001 . 

014902 (2013) 

arXiv: 1204.1409 

044906 (2014) 

arXiv:1310.8651 


arXiv:1212.1008 


[90] 

[91] 

[92] 

[93] 


22 


[94] B. Schenke, S. Schlichting and R. Venugopalan, Phys. Lett. B 747, 76 (2015) arXiv: 1502.01331 [hep-ph]]. 

[95] E. H. Lieb, “The stability of matter,” Rev. Mod. Phys. 48, 553 (1976). 

[96] G. Torrieri, Phys. Rev. C 89, no. 2, 024908 (2014) arXiv:1310.3529 [nucl-th]]. 

[97] G. Basar and D. Teaney, Phys. Rev. C 90, no. 5, 054903 (2014) arXiv:1312.6770 [nucl-th]]. 

[98] R. S. Bhalerao, J. Y. Ollitrault, S. Pal and D. Teaney, Phys. Rev. Lett. 114, no. 15, 152301 (2015) arXiv:1410.7739 
[nucl-th]]. 

[99] A. Mazeliauskas and D. Teaney, Phys. Rev. C 91, no. 4, 044902 (2015) arXiv:1501.03138 [nucl-th]]. 

[100] B. Betz, M. Gyulassy and G. Torrieri, Phys. Rev. C 84, 024913 (2011) arXiv: 1102.5416 [nucl-th]]. 

[101] Y. Tachibana and T. Hirano, Phys. Rev. C 90, no. 2, 021902 (2014) arXiv: 1402.6469 [nucl-th]]. 

[102] R. P. G. Andrade, J. Noronha and G. S. Denicol, Phys. Rev. C 90, no. 2, 024914 (2014) arXiv: 1403.1789 [nucl-th]]. 

[103] M. Schulc and B. Tomasik, Phys. Rev. C 90, no. 6, 064910 (2014) arXiv:1409.6116 [nucl-th]]. 


